//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Classe para tirar media de conjuto de dados
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include "histograma.h"
 
using namespace std;

Histograma operator* (Histograma a, double b){return a*=b;}

Histograma::Histograma()
{
	
}

//carrega um Histograma pronto direto de dados
Histograma::Histograma(const vector<vector<double> > &dados)
{
  if((dados.size() != 2) || (dados[0].size() != dados[1].size())) {
    cout << "ERRO: Input de formato errado!\n";
    exit(0);
  }

  type = 1;

  //acho o menor e maior valor
  typedef vector<double>::const_iterator it;
  hist.resize(2);
  max(dados[0][0]);
  min(dados[0][0]);
  it iter1 = dados[1].begin();
  for( it iter0 = dados[0].begin(); iter0 != dados[0].end(); iter0++){
    if ( max() < *iter0 ) max(*iter0);
    if ( min() > *iter0 ) min(*iter0);
    hist[0].push_back(*iter0);
    if(iter1 != dados[1].end()) hist[1].push_back(*iter1);
    else {
      cout << "ERRO: Input em tamanhos diferentes...\n";
      exit(0);
    }
    iter1++;
  }
  mm();
  s2();
}

//Histograma::Histograma(const vector<double> &dados=vector<double> (1,0),const vector<double> &format=vector<double> (3,0))
Histograma::Histograma(const vector<double> &dados,const vector<double> &format)
{
  type = 0;
  const vector<double> in (dados.begin(), dados.end());
  int size = 0;
  double lbox = 0;
  if(format[2] > 0) {
    lbox = (format[1]-format[0])/format[2];
    size = int(lbox);
    //lbox /= 2.;
  }
  //  hist = new vector<double> (size);

  vector<double> lhist;
  try{
  for(int i = 0; i < size; i++)  lhist.push_back(0);
  }
  catch(bad_alloc){
    cerr << "Memoria nao pode ser alocada\n";
    exit(0);
  }

  double abox = format[0];
  vector<double> nbox;
  for(int i = 0; i < size; i++){
    abox+=format[2];
    nbox.push_back(abox);
    for(int j = 0; j < in.size(); j++){
      if((abox-format[2]/2.) < in[j] && in[j] < (abox+format[2]/2.)) lhist[i]++;
    }
  }

  hist.push_back(nbox);
  hist.push_back(lhist);  
  for(int i = 0; i < format.size(); i++) config.push_back(format[i]);
  mm();
  s2();
}

// 22/11/2007 - Prepara um histograma para receber dados. Considera que grid possui a
// configuracao de um histograma, ou seja, a posicao central de cada intervalo.
// Util para intervalos nao constantes...
Histograma::Histograma(const vector<double> &grid)
{
  type = 1;

  typedef vector<double>::const_iterator it;
  typedef vector<double>::iterator it_atl;
  max(grid[0]);
  min(grid[0]);
  hist.resize(2);
  hist[0].resize(grid.size(),0);
  hist[1].resize(grid.size(),0);

  it_atl hitr0 = hist[0].begin();

  double step_med=0.0;

  for(it iter = grid.begin(); iter != grid.end(); iter++, hitr0++){
    if ( max() < *iter ) max(*iter);
    if ( min() > *iter ) min(*iter);

    *hitr0 = *iter;

    double step_up = (iter != hist[0].end()) ? (( *(iter+1) - *iter ) / 2.0)  : (( *iter - *(iter-1) ) / 2.0);
    double step_dn = (iter == hist[0].begin()) ? (( *(iter+1) - *iter ) / 2.0)  : (( *iter - *(iter-1) ) / 2.0);

    step_med+=(step_up+step_dn)/2.0;

  }
  step_med/=double(grid.size());

  //calculo largura media de cada intervalo
  config.resize(3);
  config[0]=min();
  config[1]=max();
  config[2]=step_med;
}

//Permite carregar um histograma existente para analisar!
//29/03/2007 - vou mudar essa funcao para fazer (x,y) para fazer uma que 
//carregue um histograma fazer uma funcao f(Histograma&) depois...

Histograma::Histograma(const vector<vector<double> > &dados,const vector<double> &format, double len = 0.0)
{
  type = 0;
  //acho o maior e menor valor 
  typedef vector<double>::const_iterator it;
  max(dados[0][0]);
  min(dados[0][0]);
  for( it iter = dados[0].begin(); iter != dados[0].end(); iter++ ){
    if ( max() < *iter ) max(*iter);
    if ( min() > *iter ) min(*iter);
  }

  double step = format[2];

  max(format[1]);
  min(format[0]);

  //numero de posicoes no histograma
  int np = static_cast<int>( floor( (max()-min())/step ) );
  //np++;

  hist.resize(2);
  hist[0].resize(np,0);
  hist[1].resize(np,0);

  //Calcula posicoes do histograma
  hist[0][0] = (min());

  for( int i = 1; i < np; i++)
    hist[0][i] = hist[0][i-1] + step;

  for( int i = 0; i < dados[0].size(); i++){
    
    int index = static_cast<int>( floor( 0.5+(dados[0][i]-min())/step ) );
	  
	  /*
    if( (( hist[0][index]+step/2.0 ) < dados[0][i]) && index<np ){
      cout << "# - Ajusntando +\n";
      index++;
    }
    if( (( hist[0][index]-step/2.0 ) > dados[0][i]) && index>0  ){
      cout << "# - Ajusntando -\n";
      index--;
    }
	   */
    if(index > np || index < 0) aborta("Problema fatal na funcao!");
    if( index == np ) index--;

	  // Calculando quanto do pixel esta dentro da posicao atual
	  double sup = ( ( dados[0][i]+len/2. > hist[0][index]+step/2. ) ? hist[0][index]+step/2. : dados[0][i]+len/2. );
	  
	  double inf = ( ( dados[0][i]-len/2. < hist[0][index]-step/2. ) ? hist[0][index]-step/2. : dados[0][i]-len/2. );
	  
	  double ffac = ( ( len > 0 ) ? (sup-inf)/len : 1.0 ); 
	  
    //refinando a posicao do pixel

    hist[1][index] += ( ffac * dados[1][i] ); // entra apenas quanto do pixel esta dentro da posicao no histograma
	
	  if (ffac < 1.0-1e-4) {

		  
		  // Testando quando entra nos pixels ao lado
		  int npix1 = index - 1;
		  
		  if(npix1 >= 0){
			  double up = step/2.;
			  double dn = step/2.;
			  double sup1 = ( ( dados[0][i]+len/2. > hist[0][npix1]+up ) ? hist[0][npix1]+up : dados[0][i]+len/2. );
			  
			  double inf1 = ( ( dados[0][i]-len/2. < hist[0][npix1]-dn ) ? hist[0][npix1]-dn : dados[0][i]-len/2. );
			  
			  double ffac1 = (sup1-inf1)/len; 
			  if (ffac1 > 0.0 && ffac1 < 1.0) hist[1][npix1]+= dados[1][i]*ffac1;

		  }
		  
		  npix1 = index + 1;

		  if(npix1 < hist[1].size()){
			  double up = step/2.;
			  double dn = step/2.;
			  double sup1 = ( ( dados[0][i]+len/2. > hist[0][npix1]+up ) ? hist[0][npix1]+up : dados[0][i]+len/2. );
			  
			  double inf1 = ( ( dados[0][i]-len/2. < hist[0][npix1]-dn ) ? hist[0][npix1]-dn : dados[0][i]-len/2. );
			  
			  double ffac1 = (sup1-inf1)/len; 

			  if (ffac1 > 0.0 && ffac1 < 1.0) hist[1][npix1]+= dados[1][i]*ffac1;

		  }

		  
		  
	  }

  }

  for(int i = 0; i < format.size(); i++) config.push_back(format[i]);
  mm();
  s2();

}


void Histograma::load_data(const vector<vector<double> > &dados, const vector<int> &pixmap, double len = 0.)
{
  if ((dados[0].size() != dados[1].size()) || (dados[0].size() != pixmap.size()))
    aborta("Tamanho de vetores diferentes... Reveja definicao");
  
  typedef vector<double>::const_iterator ditr;
  typedef vector<int>::const_iterator iitr;

  //Coloquei o valor do ponto junto (vector<vector<double> >) pra poder fazer um teste de 
  //intervalo. Assim posso testar se o pixmap esta realmente fazendo a coisa certa...

  ditr it_dd0 = dados[0].begin(); //fazer verificacao de intervalo
  ditr it_dd1 = dados[1].begin();
	
  // len eh o tamanho do pixel. com isso posso calcular se e quanto o pixel estravasa o seu lugar no histograma

  for(iitr it_map = pixmap.begin(); it_map != pixmap.end(); it_map++, it_dd0++, it_dd1++){
    if(*it_map >= 0){
      double up = step_up(*it_map);
      double dn = step_down(*it_map);
		
      if((hist[0][*it_map]-dn <= *it_dd0) && (*it_dd0 <= hist[0][*it_map]+up)){
		
		  
		  // Calculando quanto do pixel esta dentro da posicao atual
		  double sup = *it_dd0+len/2. > hist[0][*it_map]+up ? hist[0][*it_map]+up : *it_dd0+len/2.;
		  
		  double inf = *it_dd0-len/2. < hist[0][*it_map]-dn ? hist[0][*it_map]-dn : *it_dd0-len/2.;
		  
		  double ffac = (sup-inf)/len; 
		  /*
		  cout << "# - UP = " << hist[0][*it_map]+up << " " << *it_dd0+len/2. << endl
		  << "# - SUP = " << sup << endl
		  << "# - DN = " << hist[0][*it_map]-dn << " " << *it_dd0-len/2. <<endl
		  << "# - INF = " << inf << endl;
		  */
		  if(ffac < 0){
			  cout << "# - " << ffac << endl;			  
			  cout << "# - UP = " << hist[0][*it_map]+up << " " << *it_dd0+len/2. << endl
			  << "# - SUP = " << sup << endl
			  << "# - DN = " << hist[0][*it_map]-dn << " " << *it_dd0-len/2. <<endl
			  << "# - INF = " << inf << endl;
			  exit(0);
		  }
		  //cout << "# - " << ffac << endl;			    
		  double val = *it_dd1;
		  hist[1][*it_map]+= val*ffac;
		  
		  if (ffac < 0.95) {
			  // Testando quando entra nos pixels ao lado
			  int npix1 = *it_map - 1;
			  if(npix1 >= 0){
				  up = step_up(npix1);
				  dn = step_down(npix1);
				  double sup1 = *it_dd0+len/2. > hist[0][npix1]+up ? hist[0][npix1]+up : *it_dd0+len/2.;
					  
				  double inf1 = *it_dd0-len/2. < hist[0][npix1]-dn ? hist[0][npix1]-dn : *it_dd0-len/2.;
					  
				  double ffac1 = (sup1-inf1)/len; 
				  if (ffac1 > 0.0 && ffac1 < 1.0) hist[1][npix1]+= val*ffac1;
			  }
			  npix1 = *it_map + 1;
			  if(npix1 < hist[1].size()){
				  up = step_up(npix1);
				  dn = step_down(npix1);
				  double sup1 = *it_dd0+len/2. > hist[0][npix1]+up ? hist[0][npix1]+up : *it_dd0+len/2.;
				  
				  double inf1 = *it_dd0-len/2. < hist[0][npix1]-dn ? hist[0][npix1]-dn : *it_dd0-len/2.;
				  
				  double ffac1 = (sup1-inf1)/len; 
				  if (ffac1 > 0.0 && ffac1 < 1.0) hist[1][npix1]+= val*ffac1;
			  }
			  
			  
			  
		  }
		  //hist[1][*it_map]*=ffac;
	  }
      //else
	//aborta("Verificacao de intervalo falhou! Verifique mapa de pixels!");
    }
  }
}


/*
{


    // -------------


  int size = 0;
  double lbox = 0;
  if(format[2] > 0) {
    lbox = (format[1]-format[0])/format[2];
    size = int(lbox);
    lbox /= 2.;
  }
  //  hist = new vector<double> (size);

  vector<double> lhist;
  try{
  for(int i = 0; i < size; i++)  lhist.push_back(0);
  }
  catch(bad_alloc){
    cerr << "Memoria nao pode ser alocada\n";
    exit(0);
  }

  double abox = format[0];
  vector<double> nbox;
  vector<int> ftrack;

  for(int i = 0; i < size; i++){
     nbox.push_back(abox);

    for(int j = 0; j < in[0].size(); j++){
      bool fseek = true;
      for(int jj = 0; jj < ftrack.size(); jj++){
	if(j == ftrack[jj]) {
	  fseek = false;
	  break;
	}
      }
      if(fseek && ((abox-format[2]/2.) < in[0][j] && in[0][j] < (abox+format[2]/2.))){ 
	lhist[i]+=in[1][j];
	//cout << in[1][j] << endl;
	ftrack.push_back(j);
      }
    }
    abox+=format[2];
  }

  hist.push_back(nbox);
  hist.push_back(lhist);  
  for(int i = 0; i < format.size(); i++) config.push_back(format[i]);
  mm();
  s2();

}
*/

void Histograma::imprime(ostream& out)
{
  for(int i = 0; i < hist[0].size(); i++) out << hist[0][i] << " " << hist[1][i] << '\n';
}

void Histograma::new_data(const vector<double> &dados, const vector<double> &format=vector<double>(3,0))
{
  if(format != vector<double>(3,0)){ //ajusta se for dado vector<> format
    cout << "Ajustando ... \n";
    ajusta_config(format);
    ajusta_hist(dados);
  }

  int size = hist[1].size();
  double lbox = 0;
  if(config[2] > 0) {
    lbox = (config[1]-config[0])/config[2];
    lbox /= 2.;
  }

  double abox = config[0];
  for(int i = 0; i < size; i++){
    abox+=config[2];
    for(int j = 0; j < dados.size(); j++){
      if((abox-config[2]) < dados[j] && dados[j] < (abox+config[2])){
	hist[1][i]++;
      }
    }
  }
}

void Histograma::ajusta_config(const vector<double> &format)
{
  if(hist[0].size() == 0){ //se nao tiver nenhum dado ajusta
    for(int i = 0; i < config.size(); i++) config[i] = format[i];
  }  
  else cout << "erro: config nao ajustado\n";
}

void Histograma::ajusta_hist(const vector<double> &dados)
{
  if(hist[0].size() == 0){ //se nao tiver nenhum dado ajusta
    int size = 0;
    double lbox = 0;
    if(config[2] > 0) {
      lbox = (config[1]-config[0])/config[2];
      size = int(lbox);
      lbox /= 2.;
    }
    double abox = config[0]+config[2];
    for(int i = 0; i < size; i++){
      hist[0].push_back(abox);
      hist[1].push_back(0);
      abox+=config[2];
    }
  }
}

void Histograma::mm(void)
{
  double mm = 0;
  double ocurr = 0;
  for(int i = 0; i < hist[0].size(); i++){
    mm += hist[0][i]*hist[1][i];
    ocurr += hist[1][i];
  }

  if(ocurr > 0) m=(mm/ocurr);
  else m=mm;
}

void Histograma::s2(void)
{
  double s2 = 0;
  double ocurr = 0;

  for(int i = 0; i < hist[0].size(); i++){
    s2 += sqrt((m-hist[0][i]*hist[1][i])*(m-hist[0][i]*hist[1][i]));
    ocurr += hist[1][i];
  }

  if(ocurr > 0) ss2 = s2/ocurr;
  else ss2 = s2;
} 

inline double Histograma::media(bool atualiza=false)
{
  if(atualiza) mm();
  return m;
}

inline double Histograma::sig2(bool atualiza=false)
{
  if(atualiza) s2();
  return ss2;
}

ostream& operator<< (ostream& oo, Histograma in)
{
  
  for(int i = 0; i < in.size(); i++) oo << in.le1(i) << " " << in.le2(i) << '\n';

  return oo;
}
